2. Meshing from Neper geometries#
Neper is a free / open source software package for polycrystal generation and meshing, see https://neper.info
To run this notebook you first have to install neper.
from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
import subprocess, random
n = 100
subprocess.run(["neper", "-T", "-n", str(n), "-format", "tess"])
======================== N e p e r =======================
Info : A software package for polycrystal generation and meshing.
Info : Version 4.9.1-6
Info : Built with: gsl|muparser|opengjk|openmp|nlopt|libscotch (full)
Info : Running on 8 threads.
Info : <https://neper.info>
Info : Copyright (C) 2003-2022, and GNU GPL'd, by Romain Quey.
Info : No initialization file found (`/Users/joachim/.neperrc').
Info : ---------------------------------------------------------------
Info : MODULE -T loaded with arguments:
Info : [ini file] (none)
Info : [com line] -n 100 -format tess
Info : ---------------------------------------------------------------
Info : Reading input data...
Info : Creating domain...
Info : Creating tessellation...
Info : - Setting seeds...
Info : - Running tessellation...
Info : Generating crystal orientations...
Info : Writing results...
Info : [o] Writing file `n100-id1.tess'...
Info : [o] Wrote file `n100-id1.tess'.
Info : Elapsed time: 0.017 secs.
========================================================================
CompletedProcess(args=['neper', '-T', '-n', '100', '-format', 'tess'], returncode=0)
verts = {}
edges = {}
faces = {}
solids = []
f = open("n" + str(n) + "-id1.tess")
def randomcol():
r = random.uniform(0, 1)
g = random.uniform(0, 1)
b = random.uniform(0, 1)
return (r,g,b)
while True:
line = f.readline()
if not line: break
if line.split()[0] == "**vertex":
num = int(f.readline())
print ("found", num, "vertices")
for i in range(num):
line = f.readline()
nr,x,y,z,hhh = line.split()
verts[int(nr)] = Vertex(Pnt(float(x), float(y), float(z)))
if line.split()[0] == "**edge":
num = int(f.readline())
print ("found", num, "edges")
for i in range(num):
line = f.readline()
nr,i1,i2,hhh = line.split()
edge = Edge(verts[int(i1)], verts[int(i2)])
edges[int(nr)] = edge
edges[-int(nr)] = edge.Reversed()
if line.split()[0] == "**face":
num = int(f.readline())
print ("found", num, "faces")
for i in range(num):
l1 = f.readline()
l2 = f.readline()
l3 = f.readline()
l4 = f.readline()
nr,*_ = l1.split()
nume,*fedges = l2.split()
face = Face(Wire( [edges[int(enr)] for enr in fedges] ))
faces[int(nr)] = face
faces[-int(nr)] = face.Reversed()
if line.split()[0] == "**polyhedron":
num = int(f.readline())
print ("found", num, "polyhedra")
for i in range(num):
nr,nrf,*polyfaces = f.readline().split()
solids.append (Solid(Glue( [faces[int(fnr)] for fnr in polyfaces ] )))
solids[-1].faces.col = randomcol()
shape = Glue(solids)
geo = OCCGeometry(shape)
Draw (shape);
found 575 vertices
found 1146 edges
found 672 faces
found 100 polyhedra
with TaskManager(pajetrace=10**9):
mesh = geo.GenerateMesh(maxh=0.1)
print ("ne =", mesh.ne)
ne = 60279
Draw (mesh);